Phase diagram of the one dimensional Hubbard-Holstein Model at 1/2 and 1/4 filling 
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The Hubbard-Holstein model is one of the simplest to incorporate both electron-electron and 
electron-phonon interactions. In one dimension at half filling the Holstein electron-phonon coupling 
promotes onsite pairs of electrons and a Peierls charge density wave while the Hubbard onsite 
Coulomb repulsion U promotes antiferromagnetic correlations and a Mott insulating state. Recent 
numerical studies have found a possible third intermediate phase between Peierls and Mott states. 
From direct calculations of charge and spin susceptibilities, we show that (i) As the electron-phonon 
coupling is increased, first a spin gap opens, followed by the Peierls transition. Between these 
two transitions the metallic intermediate phase has a spin gap, no charge gap, and properties 
similar to the negative— U Hubbard model, (ii) The transitions between Mott/intermediate and 
intermediate/Peierls states are of the Kosterlitz-Thouless form, (iii) For larger U the two transitions 
merge at a tritical point into a single first order Mott/Peierls transition. In addition we show that 
an intermediate phase also occurs in the quarter-filled model. 

PACS numbers: 71.10.Fd, 71.30.+h, 71.45.Lr 



I. INTRODUCTION 

In crystalline materials where one or more of the build- 
ing blocks of the crystal structure is a large molecule the 
vibrational properties of the molecules often have large 
effects on the overall electronic properties of the mate- 
rial. One large family of such molecular crystalline ma- 
terials are the organic conductors and superconductors^. 
While some molecular crystals such as the fullerene 
superconductors^ have a three-dimensional crystal struc- 
ture, many other examples are either quasi-one or quasi- 
two dimensional, i.e. charge transport is restricted in 
certain directions due to anisotropic crystal structure. In 
addition to strong electron-phonon (e-ph) coupling to the 
molecular vibrations, electron-electron (e-e) interactions 
are often important in low dimensional materials. In this 
paper we present numerical calculations of the phase dia- 
gram for one of the simplest possible many-body models 
incorporating both these effects, the Hubbard-Holstein 
model (HHM) in one dimension (ID). In the HHM, in- 
ternal (intra-molecular) molecular vibrations are coupled 
to the local charge density of the electrons^. The elec- 
trons further interact with other electrons with an onsite 
Coulomb repulsion when two electrons occupy the same 
orbital^. Surprisingly complex effects result from this 
simple model due to the presence of both e-e and e-ph 
interactions. 

The ID HHM Hamiltonian we consider is 



h.c. 



3,v 



11 3," 



3 



l 3A n 3,l 



(1) 



where ct g (cj. CT ) are creation (annihilation) operators for 
electrons on site j with spin a, at (<Zj) are bosonic cre- 
ation (annihilation) operators for phonons at site j, and 



the electron number operator n^ a = c^ a c^ a . U is the 
Hubbard onsite e-e interaction energy u> is the disper- 
sionless phonon frequency and g is the e-ph coupling 
constant. All energies in this paper will be given in units 
of t, the electron hopping integral. 

We will concentrate primarily on EqQ]in the half-filled 
band limit (one electron per lattice site), but also discuss 
briefly the quarter-filled band (one electron per two lat- 
tice sites). The effect of e-ph interactions on a half-filled 
ID metal is well known: for inter-molecular phonons cor- 
responding to the relative motion of adjacent molecules 
in the crystal, the ID lattice dimerizes with alternating 
strong and weak bonds. In this bond-order wave (BOW) 
state the expectation value of the electron hopping be- 
tween adjacent sites alternates between strong and weak 
values. The dimerized chain then has a gap at the Fermi 
level and an insulating ground state 5 . This Peierls state 
has both charge and spin gaps, and a bond modulation 
at 2kp (q= n) at half filling. For Holstein- type phonons 
that couple to the local charge density a similar Peierls 
state occurs, but instead of bond deformation the local 
charge density is modulated in a charge density wave 
(CDW) ground state. The CDW Peierls state at half 
filling has alternating large and small charge densities 
again with periodicity 2kp. Similarly the effect of the 
Hubbard onsite interaction in ID is well known: for any 
U > at half filling, the ground state is an insulator 6 . 
Anti-ferromagnetic (AFM) spin correlations are present 
in this Mott insulating state, although no long-range an- 
tiferromagnetic order is possible in ID. At half filling, the 
2kp CDW cannot coexist with 2kp AFM correlations and 
hence the Peierls and AFM states are competing. 

Numerous previous studies have examined HHM 
within various approximations and analytic or numeri- 
cal techniques. In the limit to — > oo one can integrate 
out the phonons leaving an effective U composed of the 
sum of the Hubbard U and the effective phonon inter- 
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action, U g = U — 2g 2 /ui. For U c g > one expects the 
Mott state, while for C/ cff < one expects the Peierls 
stated. If the phonons are treated in the classical (adi- 
abatic lu —> 0) limit, one expects Peierls order for any 
g > at U = 0. However, it was shown in the spin- 
less model (Eq. [1] with a single species of fcrmion) that 
quantum fluctuations of the phonon field lead to a fi- 
nite e-ph coupling g c before the Peierls state is formed 
at U = The model with spin (Eq. [I]) has since 

been shown to also require a finite e-ph coupling for the 
Peierls transitio n 10 ' 11 . 

In addition to studies of the ID model several recent 
studies have been performed on the Holstein and HHM in 
the limit of infinite dimensions (d — oo) using dynamical 
mean-field theory (DMFT) and related method s 12 ' 13 ' 14 . 
In the d = oo model the system is metallic at U = 
also for g less than a finite value. However, an important 
distinction between d — oo results and those presented 
here is that at d — oo, the Mott insulating transition 
occurs at a finite value of U, U > 6t, while at d — 1 
it occurs for U > 0. Some similarities are found with 
our results, in particular that there is a deviation in the 
critical coupling for the Peierls transition from U e g = 
at small U—. 

Given that in the half-filled ID HHM at U = the 
ground state is metallic (no charge gap and a finite Drude 
weight) for a finite value of g, it was proposed that this 
metallic phase continues to exist between the Peierls and 
Mott insulating phases for U > O 1 ^. Subsequent nu- 
merical calculations confirmed that a metallic phase ex- 
ists for both U = and finite U— . In this paper, we 
present more detailed numerical results and analysis of 
the phase diagram. We confirm the intermediate phase 
using a different and more direct order parameter, and 
present more detailed finite-size scaling of the quantum 
phase transitions. From the finite-size dependence, we 
determine that the two transitions (Mott/intermediate 
and intermediate/Peierls) are of the Kosterlitz-Thouless 
(KT) type. We find that for larger U, the two transi- 
tions merge into a single first-order Mott/Peierls transi- 
tion. In our revised analysis, we find that the apparent 
presence of the Luttinger Liquid (LL) exponent K p > 1— 
does not imply dominant superconducting pairing corre- 
lations, but is more likely a finite-size effect. We present 
the phase diagram for three different phonon frequencies. 
We further show that at quarter filling a similar interme- 
diate phase occurs. 

The outline of the paper follows. We first give some 
details of the numerical method we used. Turning to our 
results, we discuss the U = case and then move on to 
finite U and the quarter-filled band. Finally, we conclude 
with a discussion of our data and their relation to other 
theoretical results, as well as unanswered questions for 
further study. 



II. METHOD 

We use the Stochastic Series Expansion (SSE) quan- 
tum Monte Carlo (QMC) methoS^^LSL. SSE pro- 
vides statistically exact results (no Trotter discretization 
of imaginary time is used) and has been adapted for many 
different quantum lattice models. Although this method 
has been described in detail elsewhere, we briefly describe 
here our treatment of the Holstein phonon interaction. 

In SSE, the partition function Z = Tr{e~P H } is ex- 
panded in terms of a series of sequences Sl of operators 

a S L ' »=1 

In Eq. [21 n is the length (number of operators) of each 
sequence, L the maximum allowed sequence length, and 
[3 is the inverse temperature and \a) is a basis state, here 
a direct product of electron and phonon configurations. 
In order to obtain the ground-state phase diagram, all 
results presented here used (3/t > 2N, where N is the 
number of lattice sites. The operators H ai ^ i define the 
Hamiltonian, and have type (a,) and bond (bt) indices 
with i indicating their position within the sequence Sl- 
For the ID Hubbard model (Eq. [1] with g = lu = 0), we 
have three different operators representing the diagonal 
interaction and electron hopping for both spins^: 

H W = C- -[(n TJ - ^){n u - i) 

+ ("T.j + i - \)( n U+i - 5)] 

+ /i(2 - nj - rij+i) (3) 

H2.3 = c i+i,t c J>t + h.c. (4) 

H 3,j = c ]+i,i c j,l + h - c - ( 5 ) 

Here j labels the first site of the bond the operator acts 
on. fi is the chemical potential, written here so that 
fi = corresponds to half filling. C is a constant chosen 
so that the expectation value of is always positive 
definite. In addition to the operators of Equations [5] — 
O a null operator Ho is used as a place-holder in the 
sequence expansion. We represent the phonons in the 
phonon-number basis and add the following operators for 
the e-ph interactions and phonon diagonal energy: 



Hfa = ga]n 3 (6) 

Hl 3 = ga] +1 n ]+1 (7) 

Htj = gajnj (8) 

H b.j = mj+irij+i (9) 

flflj = w{N p -a\ aj ) (10) 



Additionally, for the HHM, /1 in Eq. [3] should be replaced 
by (2g 2 /w+/i). Since the Holstein interaction couples the 
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electron density on a single site while the SSE operators 
typically act on bonds composed of two sites, we define 
two different phonon operators acting on phonon num- 
bers on the left or right of the bond. These have super- 
scripts "L" and "R" respectively. The diagonal operator 
Hq j also acts on a single site j. N p is a cutoff in the 
maximum number of phonons per site. We discuss fur- 
ther below the choice of this cutoff, but in practice it can 
be chosen large enough so as to not affect the accuracy 
of the method. 

The Monte Carlo updating is composed of an update 
for the electrons followed by an update for the phonons. 
The electron update consists of an update changing the 
number of diagonal H\j operators in the sequence, fol- 
lowed by a loop update that exchanges diagonal and off- 
diagonal operators. For the electrons we use the directed 
loop algorithm 2 ^. We note that the operators Eq. [6] 
through Eq. [TU] are not changed during the electron loop 
update. The phonon updating also consists of two parts, 
first a diagonal update changing numbers of Hqj opera- 
tors, and second an off-diagonal update exchanging Hi, 
H 4 , and H§ operators. In the diagonal phonon update, 
Hq operators are interchanged with Hgj operators with 
the following Metropolis algorithm probabilities (Nh is 
the total number of non-ifo operators present in the se- 
quence) : 



Nf3w(N p - (o]q 3 -)) 

L — Nh 
L-N H + 1 

Nf3u;(N p - (a] aj )) 



(11) 
(12) 



The phonon update for off-diagonal operators is simi- 
lar to the technique described in Reference For each 
site in the system, a subsequence is constructed which is 
a subset of the operators in Sl- The subsequence con- 
sists of only the operators Hi^ m , H 4>m , and i?5 jm which 
act on phonons at a particular site m. Within the subse- 
quence, adjacent pairs of operators are then selected at 
random and changed with a Metropolis probability. The 
pair substitutions that change the phonon number are 
(omitting the site index m as all apply to the same site) : 



(H 4 ,H 5 ) 
(H 5 ,H A ) 



(H 4 , H 5 ), (H 5 , 
{Hi, Hi) 
(Hi, Hi) 



Hi) 



(13) 
(14) 
(15) 



In addition, pair substitutions are attempted that swap 
the order in the subsequence of the two operators. When 
two different pairs may be substituted the substitution 
made is chosen randomly. Note that the L and R indices 
in Eq. [6] through Eq. [9] are not needed during the pair 
updating, but updates involving the Hi operators must 
be canceled with 50% probability (for each Hi operator 
in the pair). If a Hi operator changes into a phonon 
operator as a result of the update, a L or R index is 
assigned when the subsequence update is completed and 
merged into Sl- The Metropolis substitution probabili- 
ties depend on phonon as well as diagonal electron matrix 



elements variables, the e-ph coupling constant g, and the 
number N<i of diagonal phonon operators (Hqj) that are 
present between the two operators of the pair. N^ may 
be stored when the subsequence is constructed. For ex- 
ample, in terms of just the change in the phonon part of 
the operator, 



P[(Hi,Hi) -» (H 5 ,H 4 )} = Rng 2 



P[(H 5 ,H 4 )^(Hi,Hi)} = 

ng z 



N p -n + l 

N p -n 

N p -n 
N p -n + l 



(16) 



N,i 



(17) 



where n is the number of phonons present in the sequence 
position just before the operator pair. R in Eq. [16] and 
Eq. [T7] is the ratio of diagonal matrix elements from the 
electronic Hamiltonian. In practice, the number of pair 
substitutions performed is chosen to be approximately 
the same as the number of operators in the subsequence. 

We use standard methods to calculate various observ- 
ables within our SSE code^i. To determine phase bound- 
aries of the model we primarily use the charge and spin 
susceptibilities at wavevector q given by 



(18) 



XpAl) = 



N ^ 

3,k 



"^- fe > / dr(OHr)Otm (19) 



In Eq. [TUthe charge susceptibility Xp(l) (spin suscepti- 
bility X(r{q)) corresponds to the + (— ) sign. Similarly, 
we also use the static structure factors, S p (q) and S a (q): 



1 

N 



3,k 



- k) (ofo± 



(20) 



The effective low-energy properties of many interacting 
ID models can be understood in terms of a LL picture, 
and the asymptotic properties of the system described 
by a small number of parameter o 23 i 24 . In particular, the 
asymptotic decay of correlation functions can be related 
to the correlation exponents K p and K a for charge and 
spin respectively. In the long wavelength limit, these ex- 
ponents may be calculated from the slope of the structure 
factors: 



Kn 



7r lim Sp Jq)lq 

q— »0 



(21) 



In practice, one uses the behavior of nS(q) j q at the small- 
est available q for the periodic ring, q\ = 2n/N . With 
proper finite-size scaling in N , this gives the Luttinger 
liquid exponent for the system 2 ^. Based on calculations 
of acoustic phonons coupled to ID electrons it has been 
suggested that the expected relationship of K p to the 
correlation functions must be modified in the presence 
of phonon interactions with retardation 2 ^ 2 ' 28 . We will 
discuss this further in Section fill Dl below. However, we 
note that the interpretation of K a is not modified in the 
presence of phonon retardation effects since spin-rotation 
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FIG. 1: (color online) Slope of the spin structure factor at 
wavevector qi = 2n/N versus e-ph coupling g for a 16 site 
half- filled system with (7 = 4 and uj — 1. nS a (qi) /qi crossing 
one indicates the opening a spin gap. Different symbols show 
the convergence with increasing phonon cutoff N p . 



symmetry is preserved in the HHM. K a is expected to be 
exactly equal to one unless a spin gap is present, and the 
condition that TtS a {q\) / 'q\ decreases below one is a sensi- 
tive indicator for the opening of a spin gap (sec Fig.[l} 20 . 
We find that finite-size effects in determining the phase 
boundaries using Eq. [21] are worse than when using the 
susceptibilities, Eq.[19j due to the necessity of taking the 
limit qi — > in Eq. 1211 Therefore we will primarily use 
the susceptibilities in order to determine the phase dia- 
gram boundaries. 

We choose the phonon cutoff -/V p such that phonon 
occupation numbers during the simulation never reach 
within some fraction (~20%) of the cutoff, similar to the 
method in which the maximum sequence length L is set 
self-consistently in SSE simulations. We have verified 
that our results are converged with respect to N p . Typ- 
ical variation with iV p is shown in Fig. [T] for a 16 site 
system with U = 4 and u) = 1. We find that choosing 
N p too small can have a noticeable effect on the criti- 
cal coupling for transitions, and especially on quantities 
measured in the Peierls phase. 

Autocorrelation time r is an important measure of 
the overall efficiency of a Monte Carlo method. Cor- 
relations between measurements are typically found to 
decay as ~ e~*/ r , where t is in units of Monte Carlo 
time corresponding to the number of updating steps com- 
pleted. If measurements are correlated the estimated sta- 
tistical error must be increased. In general, it is found 
that near quantum phase transitions, r often increases 
steeply, making calculations near phase boundaries diffi- 
cult or impossible. One tool available to improve QMC 
calculations near phase boundaries is quantum paral- 
lel tempering 29 . In this technique, separate processors 
on a parallel computer have slightly different parame- 
ter values. Periodically a Metropolis move is attempted 
to switch the configuration between adjacent processors. 
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FIG. 2: (color online) Integrated autocorrelation time for 
N = 16, U = 2, u> = 1, f3 = 32, N p = 30, as a function of e-ph 
coupling. Filled circles (squares) are the autocorrelation time 
for the charge (spin) structure factor S(qi) at qi — 2-k/N. 
Open symbols are for the same observables, but calculated 
using quantum parallel tempering. Arrows indicate the lo- 
cation of the two transitions (see Section ITTT)) . We find that 
parallel tempering significantly reduces the autocorrelation 
time near the transitions. 

These moves help to prevent the algorithm from getting 
"stuck" in one configuration, and consequently reduce 
the autocorrelation time. In Fig. [2] we show the inte- 
grated autocorrelation time for long wavelength structure 
factor measurements (Eq. [2~T1) . defined as in Reference l2~il. 
Our definition of one Monte Carlo step is similar to Ref- 
erence Hi], with an average 2Nh loop vertices visited in 
the electronic loop update. As expected we find that 
t increases greatly near the Peierls transition. We also 
find that parallel tempering decreases the autocorrela- 
tion time significantly and is essential to obtain reliable 
results near the Peierls transition. 



III. RESULTS AT HALF-FILLING 

We first present our results for the half-filled band, first 
in the case U — and then for finite U. 



A. U — : the Peierls transition 

Eq. Q] has been studied in great detail for the case of 
U = 0. One of the key questions is whether the transition 
to the Peierls state occurs for finite critical coupling or 
for any value of g > 0. The transition occurring at finite 
g is expected to be of the KT typei^. KT transitions 
at finite phonon coupling have been found in a number 
of ID phonon-coupled models including the spinless Hol- 
stein model (Eq. [1] with only one species of fermion)^^, 
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FIG. 3: (color online) Finite-size scaling of the q — n charge 
susceptibility for (7 = and uo = 1 at half filling. Data are 
for system sizes up to N = 128 sites. In (a) we plot Xpi^/N 
versus N. At critical coupling Xp( n )/N approaches a constant 
for large N. Note that the g = curve corresponds to free 
fermions (no phonons). In (b) Xpi^/N is plotted versus the 
effective e-ph coupling 2g/u, for system sizes N — 8 (open 
circles), 16, 32, 64, and 128. The inset shows finite-size scaling 
of the transition point obtained by plotting the value of 2g 2 /o; 
where Xp( n )/N for system size N exceeds the susceptibility 
for system size N/2. Line in inset is a linear fit. We estimate 

1.00 (g c2 « 0.71). 



the critical coupling as 2gi 2 /u 



the XY model coupled to dispersionless phonons^, the 
Heisenberg model coupled to dispersionless phonons 19 , 
and the extended Peierls-Hubbard model coupled to dis- 
persionless bond phonons- . We confirm that indeed a 
finite critical coupling exists and show that the finite- 
size scaling of the observables is consistent with a KT 
transition. 

A KT quantum phase transition is difficult to de- 
tect because the gap opens exponentially slowly. For 
Holstein-type phonons that couple to the local electron 
density, the appropriate order parameter for the transi- 
tion is the 2kp charge susceptibility. The critical coupling 
(we will denote the critical g for the Peierls transition 
as g C 2) may be determined from the finite-size scaling of 




FIG. 4: (color online) Long-wavelength spin and charge struc- 
ture factor slopes for U — and u) = 1 at half filling. Open 
(filled) symbols are for charge (spin). Data are for system sizes 
of N = 16, 32, 64, and 128 sites. For any g > TrS a (qi)/qi 
is less than 1 indicating the presence of a spin gap. The in- 
set shows the finite-size scaling of the point where S p (qi)/qi 
crosses one, with the line a fit to a quadratic. We estimate 
the critical coupling as 2g\ 2 /io ps 0.85. The appearance of 
these data for g < g c2 are similar to those for the negative 
U Hubbard model (Fig. [7] for U < 0). The interpretation of 
these data is discussed in Section MIDI 



the 2kp charge susceptibility, Xpi n )- Xp( n )/N should ap- 
proach zero logarithmically below g c2 and should diverge 
above g C 2- Exactly at g = g C 2, log corrections vanish and 
X P {ir)/N should approach a constant value with increas- 
ing N. Our SSE results confirm that x p (7r) does scale 
in this manner. In Fig. |3)Ja) we show charge susceptibil- 
ity data for U = and oj = 1, which is consistent with 
a KT transition at g C 2 ~ 0.7. We see a clear decrease 
of Xp{^)/N with system size below the transition and a 
clear increase above the transition. Plotted as a function 
of effective e-ph coupling 2g 2 fw (Fig. |3Jb)), XpW/N for 
different TV cross at the transition. In the Fig. [UJb) we 
show a finite-size scaling of the transition point obtained 
by plotting value of 2g 2 /u where the susceptibility curve 
for N sites intersects the data for N/2 sites. We find 
that these intersection points are well fit to a linear de- 
pendence in 1/N, giving 2flr^ 2 /a»=1.00 for U = 0. 

In Fig. U we show for comparison the long- wavelength 
charge and spin structure factor slopes, Eq. [2B which 
are estimates for the LL exponents K p and K a . For any 
g > K a is less than 1 and decreases with increasing 
chain length, indicating a spin gap. Furthermore, in the 
spin susceptibility (not shown here), we find no sign of 
any transition at the critical coupling where x P (7r)/7V 
diverges. We denote the critical coupling for the spin 
gap opening as g c \ . Hence we conclude that a spin gap is 
present for any g > g c \ = when U = 0, but a charge gap 
is only present for g > g C 2- In the inset of Fig. [4] we show 
the finite-size scaling of point where K p — 1 (2g 2 2 /uj). 
We discuss further the K p data in Section IIIIDI and the 
apparent small discrepancy between g C 2 determined from 
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FIG. 5: (color online) (a) Charge (open symbols) and spin 
(filled symbols) susceptibility for the half-filled HHM with 
U = 2, uj = 1. The first transition (g c i) occurs where 
Xp — Xo-, corresponding to U e g = 0. The second transition 
(gc2) is the Peierls transition, where Xp( n )/N diverges as in 
Fig. [3] Note that the spin susceptibility is also divided by N 
to make the crossing at [7 c ff = clear, (b) Long-wavelength 
spin structure factor for U = 2, w = 1. The point where 
TrScr(qi)/qi crosses unity indicates the opening of the spin 
gap, identical to the point where \p = X" in ( a )- 

susceptibility versus K p data. 

B. U > : Intermediate phase 

We next consider the case with U > at half filling. 
To avoid any possible difficulties of interpreting numeri- 
cal estimates for K p , we determine all phase boundaries 
directly from susceptibilities and K a . In the ID Hub- 
bard model (g — in Eq. [l} charge and spin degrees of 
freedom effectively switch places at U = 0. In terms of 
the susceptibilities, Xpi^) an d Xo-M are exactly equal at 
U = 0. 

In Fig. [5fa) we first show the 2kp charge and spin 
susceptibilities for U = 2 and u> = 1. We find that when 
[/ c ff w (g = 1 for U = 2 and oj = 1), the charge and spin 
susceptibilities become equal as in the simple Hubbard 
case. The estimate for K a , shown in Fig. EJb), again 
crosses one indicating an opening of a spin gap. This 
transition is therefore the same transition g c i as discussed 
above in Section TlII Al but now occurring at finite g. The 
quantum phase transition as g increases past g c \ appears 
identical to the transition as U becomes negative in the 
ID Hubbard model. Based on the similarity with the ID 
Hubbard model, we conclude that the spin gap transition 
here is also of the KT form. 



In Fig. EJa) a second transition takes place beyond 
the spin gap transition at g c \. This second transition is 
again the Peierls transition indicated by the divergence of 
X p (tt)/N. Beyond the second transition point (<? > g C 2) 
X p (tt)/N increases with increasing system size, and as 
in Fig. [3fb) Xp( 7r )/-^ f° r different system sizes cross at 
.9 = <7c2 when plotted versus e-ph coupling. For g c i < 
g < 9c2 we now have a third intermediate phase, which 
has a spin gap but no Peierls order. In Fig. [5la) we see 
only very small finite-size effects in determining g c \ and 
g C 2 from the susceptibility data. The g c \ from our data 
shows little deviation from U e s = 0, at least for small 
to intermediate U as compared to u>. Finite-size effects 
are more significant in K a as estimated from the spin 
structure factor slope in Fig. [5]Jb) because q\ = 2tt/N 
only approaches q± = in the limit N — > oo. However, 
for increasing N, g c \ as estimated from K a does converge 
to the same value we obtain from the susceptibility. 

As U increases, we find that two transitions at g c \ 
and g C 2 occur closer together, becoming indistinguishable 
from each other at approximately U ~ 5 for uj = 1. At 
this point and for larger U, the two KT transitions merge 
into a single Mott/Peierls transition. We next show that 
this merged transition is first order. 



C. First order transition 

Above a critical U value U — U m we find that the spin- 
gap and Peierls transitions coincide. The phase diagram 
then has a shape very similar to that of the half-filled 
ID extended Hubbard model (EHM)^^^^. In the 
half-filled EHM, as the nearest-neighbor Coulomb repul- 
sion V is increased for fixed U, there is a transition from 
AFM to CDW order. This transition is continuous for 
small U and first-order for U > U m . In a first-order 
quantum phase transition, observables become discon- 
tinuous as one of the Hamiltonian parameters is varied. 
For the HHM, a change to first order behavior for strong 
coupling has also been seen in DMFT studies^. As m 
Referenced we take [m(p)] 2 = S p (n)/N as an order pa- 
rameter for the Peierls CDW state. We show in Fig. [5] 
[m(p)] 2 for U — 8 and oj = 1. We find a sharp jump in 
[m(p)] 2 at the transition with the discontinuity becoming 
stronger for larger system sizes. Other observables such 
as the ground state energy and bond order also show 
discontinuous behavior consistent with a first order tran- 
sition. In fact, this point is a multi- critical point. In 
the EHM there is an intervening phase with long-range 
BOW for U < Urrir2<&2!L. We find very similar behavior 
in the HHM except that the intervening phase here is the 
metallic intermediate state. We cannot calculate a pre- 
cise value for U m , but for u) = 1 it appears comparable 
(U m ~ 5 for u> = 1) to the value found in the half-filled 
EHM, U m = 4.7±0.l2& We also remark that the change 
in the order of the transition may be related to discus- 
sions of quantum to classical crossover in e-ph coupled 
models^. 
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FIG. 6: (color online) First order Mott/Peierls transition for 
large U. We show the CDW order parameter [m(p)} 2 
text) vs. —U e ff = 2g 2 /uj — U for oj = 1 and [7 = 8. 
Mott/Peierls transition occurs for U e g ~ —0.2. 
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FIG. 7: (color online) LL exponents for the ID Hubbard 
model (Eq.[TJwith g — 0) estimated from the long-wavelength 
charge and spin correlations. K p (K a ) is given by open (filled) 
symbols. In the infinite N limit, K p — 1 for any U < and 
K p = for U > 0; K a = 1 for any 17 > and *f CT = for 
any [7 < 0. Observing these limiting values (shown by full 
and dashed horizontal lines) is difficult due to finite value of 
qi — 2ir/N and also the logarithmic scaling with N for the 
exponent whose value is unity. 



D. Discussion of Luttinger exponents 

In the LL picture, K p and K a determine the asymp- 
totic decay of correlations functions and hence measure- 
ments of these exponents in finite systems have often 
been used to determine the phase diagrams of ID models. 
Specifically, K p > 1 corresponds to attractive charge cor- 



relations, while K p < 1 corresponds to repulsive charge 
correlations. It is first instructive to review the LL ex- 
ponents for the ID Hubbard model and sources of error 
in finite-size systems. At half filling for U > the ID 
Hubbard model is insulating (K p — 0) with no spin gap 
(K a = 1, spin rotational invariance holds). For U < 
there is a spin gap (K a = 0), and degenerate CDW and 
singlet superconducting (SS) pair correlations [K p = 1). 
Therefore, the LL exponents are discontinuous at U = 0. 
The transition at U = is of the KT type, with the gaps 
(charge gap U > or spin gap U < 0) opening exponen- 
tially slowly as U is varied from zero. In Fig. [7] we show 
K p and K a for the ID Hubbard model calculated using 
Eg. 1211 There are two primary sources of finite-size error: 
first, the requirement that q — > in Eq. 1211 and second, 
the presence of logarithmic scaling corrections near a KT 
transition. The scaling with system size is slow close to 
the transition (U = 0) and particularly slow for the ex- 
ponent that is expected to be equal to 1 (K a for U > 
and K p for U < 0). Such logarithmic scaling has been 
noted in other ID electron and spin models and makes it 
difficult in practice to observe K a = 1 for the positive-C/ 
Hubbard model in a finite-size calculatio n 19 ' 20 . As dis- 
cussed in Section IIII Al log corrections are expected to 
vanish exactly at critical coupling. In Fig. [7] this occurs 
at U = 0, where K p and K a curves for all system sizes 



cross at K n = K„ 



1. 



Turning now to the HHM, the variation of K a for 
.9 < Sci (Fig. Hlb)) is consistent with log corrections 
in the spin degree of freedom that vanish at the spin 
gap transition. This observation further reinforces our 
statement that the spin gap transition is also of the KT 
type. For the K p data in Fig. 01 K p at g = is again 
exactly unity. K p then crosses one from above at a g 
roughly consistent with the g c i determined from the sus- 
ceptibility data in Fig. [3] Assuming the Peierls transi- 
tion occurs where K p = 1 gives a critical coupling of 



2gi 



0.85 after performing finite-size scaling using ./V 



up to 128 sites (Fig. HJ. The form of the K p plot for the 
HHM (U — 0) is clearly similar to K p for the negative 
U Hubbard model (Fig. [7J, with K p starting at one for 
zero coupling, and becoming slightly larger than one for 
nonzero coupling. While this apparent K p > 1 may be 
interpreted as meaning that superconducting pair corre- 
lations are dominant^, a more plausible interpretation 
is that the apparent K p > 1 is a consequence of logarith- 
mic scaling corrections. This implies that the true K p 
should be exactly equal to unity for g < g C 2, and drop to 
zero for g > g C 2- This further implies that the intermedi- 
ate state has degenerate CDW and SS correlations. This 
statement is consistent with our finding that the U = 
HHM for g < g C 2 has a spin gap but no charge gap. 

Calculations for a model of acoustic phonons coupled 
to ID electrons found that the LL expressions for decay of 
correlation functions must be modified due to retardation 
effects^. Specifically, the dominance of CDW and SS 
correlations is given by 



K P A 



< 1 



(CDW) 



(22) 
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B/K p < 1 (SS) 



(23) 



where A and B depend on the strength of the e-ph 
coupling 26 . With zero e-ph coupling, A = B = 1. For 
increasing e-ph coupling, A > 1 and B < 1, with A di- 
verging and B approaching a finite value. The renormal- 
ized boundary for the metallic/Peierls transition is then 
K p = 1/A. While there is no reason to expect that for 
the HHM model (with dispersionless phonons) the LL re- 
lations should be renormalized in the same manner, our 
SSE data may be consistent with 1/A slightly less than 
1. Upon close examination of Fig. [4] and Fig. [3l the g c i 
as determined by K p crossing one is slightly smaller than 
the g C 2 determined by susceptibility. The g C 2 determined 
from K p (Fig. 3]) would coincide with the g C 2 determined 
from Xp( 7r ) (Fig- Elk)) if the horizontal line in Fig. [4] is 
moved slightly below one, or 1/A ss 0.95. 

For larger U, the size of the intermediate region 
shrinks, and K p peaks at the transition, with K p ap- 
proaching one with increasing N (see Fig. 2(a) in Refer- 
ence HH) U — 2, uj = 0.5). The peak at the transition is 
consistent with K p = in the Mott and Peierls states, 
and K p = 1 only along their boundary. The apparent 
K p < 1 at the peak may be due to the closer proxim- 
ity to the first order transition, where K p drops quite 
rapidly to zero. For U — 2 and u> = 0.5, we estimate 
that 0.95 < 1/A < 1. If renormalization as in Reference 
l26l does occur, for all parameter values we investigated it 
appears that the effect is relatively small (0.9 < A < 1). 
Because measuring SS correlations is not practical in the 
SSE method, we cannot determine a value for B. Eq. [23] 
with B < 1 would imply that a SS correlations are dom- 
inant whenever K p exceeds a value that is smaller than 
one. SS is dominant for any nonzero e-ph coupling for 
U = in the calculation of Reference [2y, which seems 
unlikely in the HHM. We will discuss these implications 
further in Section M 



E. Phase diagram, half filling 

In Fig.[8]we show the phase diagram for to = 0.5, u> = 1, 
and uj = 5. All points were determined using suscepti- 
bility data for systems up to 32 (and in some cases 64 
and 128) sites. We find that with increasing u> the width 
of the intermediate region increases, and the tricritical 
point U m moves to larger U. One further observation 
is that for U > U m , the deviation of the Mott /Peierls 
boundary from U e s = becomes noticeable, with the 
boundary shifting to C/ cff < (above the dotted lines in 
Fig. [8]). This shift can be seen for example in Fig. [6l For 
U < U m , the Mott/intermediate spin gap boundary is 
very close to the line U c s = 0. 



IV. QUARTER FILLING 

Many of the materials that the HHM is most applicable 
to are not half-filled. For example, most of the quasi-lD 
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FIG. 8: Phase diagram of the half-filled HHM for u = 0.5, 
u — 1, and uj = 5. The dashed line is given by U — 2g 2 /cu. 
All phase boundaries are determined using susceptibility and 
K a data, with uncertainty approximately the size of the sym- 
bols. Lines are guides to the eye. The three phases shown are 
Mott, (I)ntermediate, and Peierls. The Mott/I and I/Peierls 
boundaries merge into a single first-order Mott/Peierls bound- 
ary indicated by a heavy line for U > 4 for uj = 0.5 and U > 5 
for uj = 1. 



organic superconductors are 3/4 filled (1/4 hole filled) 1 . 
We therefore present some results for the HHM at quarter 
filling. Although for many of these materials it is neces- 
sary to include long ranged Coulomb interactions (the 
extended Hubbard V term)^, we will continue to focus 
on the HHM Hamiltonian with only onsite U and e-ph 
terms. We comment on the expected effect of V further 
below. As quarter filling is commensurate a Peierls state 
is also expected to occur for sufficiently large g. There 
are however significant differences between half-filled and 
quarter-filled Peierls states. At quarter filling there are 
more than one possible pattern of charge and bond dis- 
tortion, and which one actually occurs depends on the 
values of U as well as U 39 i 40 . In the absence of phonons, 
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O-O 2g7(B=1.69 
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FIG. 9: (color online) (a) Charge (open symbols) and spin 
(filled symbols) susceptibilities for the quarter-filled HHM 
with U = 2, u = 0.5. (b) Long-wavelength spin structure 
factor for same parameters. We find similar behavior to half 
filling, Fig. O with first a transition to a spin-gapped state, 
and second the transition to the Peierls CDW state. 



the quarter-filled band for finite U is a LL with neither 
charge nor spin gaps. As at half filling, Xp(^f) and 
Xa(2kp) are degenerate at U = (note that 2kp = 7r/2 at 
quarter filling and corresponds to a correlation function 
with period 4 in real space). In the presence of phonons, 
we again expect the charge susceptibility Xp(^f)/N to 
diverge. 

Our SSE results show that the HHM at quarter fill- 
ing is in many respects similar to the half-filled case. 
In Fig. [5] we show Xp{^f)/N and nS cr (qi)/qi versus 
2g 2 /u. We again find two transitions: first a transi- 
tion to a spin-gapped state, and second the transition 
to the Peierls state. As at half filling, the spin gap opens 
very close to the point where U e s = U — 2g 2 /ui = 0. 
The phase diagram at quarter filling is therefore nearly 
identical to the phase diagram at half filling, with LL, 
intermediate, and Peierls phases. We find that the in- 
termediate phase is slightly wider at quarter than half 
filling. For example, at quarter filling with U = 2 and 
u = 0.5 (Fig. 0, 2g%_/b> w 1.7 and 2g 2 c2 /u w 2.6, com- 
pared to 2g 2 1 /u> w 2.0 and 2g 2 2 /Ld w 2.3 at half filling. 
We note (see Fig. [9j> that at quarter filling we see slightly 
greater deviation from the U c g = in the first (spin gap) 
transition. At present we do not have enough SSE data 
to investigate whether the tricritical point U m occurs as 
at half filling, but in our data at quarter filling we do 
find that with increasing U g c \ and g c i become closer to- 
gether. This suggests that a tricritical point also exists 
at quarter filling. 
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FIG. 10: (color online) Charge-charge correlations {rurij) ver- 
sus distance \i — j\ for a 32 site quarter-filled system with 
U = 2 and u = 0.5. The three values of 2g 2 /uj = 1.69, 2.25, 
and 2.89 correspond to LL, intermediate, and Peierls states, 
respectively. We find that in all three regions the charge cor- 
relations at quarter filling are of the form • • -2000- ■ ■. 



At quarter filling there are two possible CDW's that 
are period-4 (2kp). These have charge densities in car- 
toon form of either ■ • ■ 1 100- • • or ■ • -2000- • •, where "1" 
or "2" indicates a charge density greater than the aver- 
age density of 0.5 and "0" indicates a charge density less 
than the average^ 9 -. The pattern • • -2000- ■ • is found in 
the uncorrelated (U — 0) band. In Fig. [TO] we plot the 
real-space charge-charge correlation function (niTij) ver- 
sus distance \i — j\ for a range of </'s in the three phases. 
We find that the charge-charge correlation function peaks 
for sites separated by four lattice sites, consistent with a 
CDW state of the • ■ -2000- ■ • form. The strength of the 
CDW correlations does not greatly change going from the 
LL to the intermediate phase, but increases rapidly after 
the Peierls transition. In the • • -2000- • • CDW, the three 
small charges are not exactly equal, and the actual charge 
densities are in sequence large, medium, small, medium 
(LMSM). This charge pattern coexists with a BOW be- 
cause L-M and M-S bonds are inequivalent. Fig.fTOlshows 
that the charge correlations follow this LMSM pattern as 
expected. We conclude that the pairing at quarter filling 
in the HHM consists of onsite electron pairs as found at 
half filling, at least for the small through intermediate U 
we have currently investigated. 

The distinction between these two CDW patterns at 
quarter filling is important because while • ■ -2000- • ■ is 
related to onsite electron pairs, the more extended CDW 

• • -1100- • • is related to nearest-neighbor pairing. The 

• • -1100- • • requires bond-coupled phonons in addition to 
the Holstein phonons considered her o 39 ' 40 . In addition, 
the pattern of the BOW (the location of the "strong" 
bond) coexisting with the • • -1100- ■ ■ CDW also depends 
on the strength of V—. If a similar metallic phase exists 
adjacent to the ■ ■ ■ 1 100- • • CDW, it is possible that a re- 
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gion of nearest-neighbor superconducting pairing will be 
found that may be relevant to real quarter-filled molecu- 
lar superconductors. 

V. CONCLUSIONS 

To summarize, we have presented numerical data for 
charge and spin correlations of the ID HHM model at 
half and quarter filling. We have based our phase dia- 
gram on charge and spin susceptibilities, which provide 
direct indication of phase boundaries with much weaker 
finite-size effects than previous calculations based on LL 
exponents^ 6 -. We find that the spin gap and Peierls tran- 
sitions do not occur simultaneously unless U is larger 
than a critical U m . For U < U m as the e-ph coupling is 
increased from zero, the spin gap opens before the Peierls 
state forms. The intermediate state is metallic with a 
spin gap but no charge gap, and the transitions to and 
from the intermediate state are of the KT type. Our 
physical picture of the intermediate state is that at the 
spin gap transition (g c i), pairs are formed, but are dis- 
ordered and do not order in a Peierls state until the e-ph 
coupling is further increased. For U > U m , the two tran- 
sitions merge into a single first-order Mott/Peierls transi- 
tion. With finite-size calculations we cannot completely 
discount the possibility of a small charge gap (small com- 
pared to the finite-size gap) in the intermediate region. 
However, finite charge stiffness (Drude weight) provides 
further evidence for metallic behavior in the intermediate 
stated. 

Comparing to other calculations, the critical coupling 
we determined for g C 2 at U = is consistent with previ- 
ous result s 10 ' 11 . The variational results of Reference [ID 
find the intermediate phase existing in a narrow region 
on both sides of the U e g = line, while we find the inter- 
mediate phase only for U e g < 0. Several calculations of 
the single-particle spectral function are available for the 
HHMiiii^, the spinless Holstein model 4 ^, as well as 
the d = oo studies previously mentioned. In Reference 42 
using a cluster perturbation theory method applied to the 
ID HHM, a small nearly dispersionless peak was found in 
the spectral function for small k. This small peak is also 
found in the spectral function of the metallic phase of the 
spinless Holstein model and may possibly be associated 
with the intermediate phase. 

Considering the possible modification of the LL equa- 



tions in the presence of retarded e-ph interactions, we find 
that while this could possibly occur in a form that would 
agree with Reference^, the amount of renormalization is 
small (1/A ~ 0.95), and possibly within finite-size errors 
in our determination of the transition points. We also 
do not see any measurable or consistent change in the 
constant A when comparing u — 0.5 and u = 1, which 
would be expected to change the amount of retardation 
in the e-ph interaction. We are not able to calculate a 
value for B in Eq.[53] However, if Eq.[53]is correct for the 
HHM model with B < 1, SS correlations would actually 
be enhanced because of retardation 26 . 

An important question is the strength of SS correla- 
tions within the intermediate region. In terms of the 
LL framework, models with K p > 1 have dominant su- 
perconducting correlations. Indeed, our numerical data 
appears to show K p > 1 in the intermediate region, but 
this result is likely to be a finite-size effect. If we set aside 
any renormalization of K pi our conclusion based on com- 
parison with the ID negative-?/ Hubbard model is that in 
the intermediate region, K p is exactly equal to one. This 
implies that in the intermediate region that CDW and 
SS correlations are in fact exactly degenerate. This exact 
degeneracy may not be easily observable in a finite sys- 
tem due to the finite size difficulties near KT transitions. 
As the SSE method is based on a world-line approach in 
imaginary time, there is no simple way to measure cor- 
relations involving four particles, which would be needed 
to measure SS correlations directly. Other QMC— and 
DMRG~ calculations suggest that while SS and CDW 
correlations are nearly degenerate in the intermediate re- 
gion, CDW correlations appear slightly stronger at long 
range. DMRG results for a ladder system suggest that 
going beyond ID can break the degeneracy, giving a re- 
gion with dominant SS correlations 2 ^. Finally, we remark 
that in addition to the logarithmic scaling difficulties, ob- 
serving metallic behavior in close proximity to a CDW 
is difficult due to the typically rapidly increasing auto- 
correlation time for QMC methods. We do find in our 
method that the autocorrelation time r increases rapidly 
close to the Peierls boundary, but believe that the use of 
parallel tempering can reduce r enough to obtain reliable 
results. 
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